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ABSTRACT 

We present a new technique for choosing spatial regions for X-ray spectroscopy, called "con- 
tour binning". The method chooses regions by following contours on a smoothed image of the 
object. In addition we re-explore a simple method for adaptively smoothing X-ray images ac- 
cording to the local count rate, we term "accumulative smoothing", which is a generalisation 
of the method used by FADAPT. The algorithms are tested by applying them to a simulated 
cluster data set. We illustrate the techniques by using them on a 50 ks Chandra observation 
of the Cassiopeia A supernova remnant. Generated maps of the object showing abundances 
in eight different elements, absorbing column density, temperature, ionisation timescale and 
velocity are presented. Tests show that contour binning reproduces surface brightness con- 
siderably better than other methods. It is particularly suited to objects with detailed spatial 
structure such as supernova remnants and the cores of galaxy clusters, producing aesthetically 
pleasing results. 
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1 INTRODUCTION 

Many X-ray telescopes such as XMM-Newton, Chandra and 
ROSAT contain detectors which capture individual photons, record- 
ing their energy and position on the detector. Therefore unlike con- 
ventional optical observing techniques, we simultaneously collect 
imaging and spectroscopic data for each part of the object. In addi- 
tion the Chandra X-ray Observatory has very high spatial resolu- 
tion (~ 1 arcsec), allowing the properties of the emitting object to 
be studied in unprecedented detail. 

With this spectroscopic and imaging information we can select 
"events" from the observation corresponding to a particular part 
of an object with a "region filter". From these events a spectrum 
can be built-up. An X-ray spectral package such as XSPEC (Amaud 
1996) can then be used to fit a physical model to the spectrum. 
Conventionally simple geometric shapes such as annuli, sectors, 
boxes or ellipses are used to define a region filter. Using sectors, 
for example, and assuming spherical or elliptical symmetry, one 
can account for projection in a cluster of galaxies. However most 
extended objects are not symmetric when observed in detail (for 
example, the Perseus cluster [Fabian et al 2000; Sanders et al 2004], 
the Cassiopeia- A supernova remnant [Hwang et al 2004], and Abell 
2052 [Blanton, Sarazin & McNamara 2003]). 

Given the morphological diversity of extended X-ray sources 
it is important to have techniques which allow us to analyse the 
spectral variation of a source over its extent. We first investigated 



E-mail: jss@ast.cam.ac.uk 



this problem when looking for cool gas in a sample of X-ray clus- 
ters using ROSAT PSPC archival data (Sanders, Fabian and Allen 
2000). We devised a technique which used adaptively smoothed 
maps (Ebeling et al, in preparation) to define contours in surface 
brightness. The ratio of the number of counts in different energy 
bands between each contour was used to define an X-ray colour. 
By using a grid of models, the absorption and temperature of the 
gas could be estimated between the contours. 

We approached this problem again with the advent of data 
from Chandra with its high spatial resolution. We created an al- 
gorithm called "adaptive binning" (Sanders & Fabian 2001) which 
used the uncertainty on the number of counts or the error on the 
ratio of counts in different bands to define the size of binning re- 
gion used. The process was simple: pass over the image, copying 
those pixels which have a small enough uncertainty on the number 
of counts or colour to an output image. Bin up the remaining pix- 
els by a factor of two. Repeat until all the pixels have been binned. 
On the final pass we bin any pixels which are not yet binned. This 
simple approach works well and was used by us and other authors 
on data from several clusters of galaxies (e.g. Centaurus - Sanders 
& Fabian 2002; Perseus - Fabian et al 2000; Abell 4059 — Choi et 
al 2004). 

The disadvantage of this approach is that the binning scale 
varies by factors of two. It is very noticeable where the scale 
changes, and some regions are overbinned. Therefore we started 
using the "bin accretion" algorithm of Cappellari & Copin (2003). 
The algorithm adds pixels to a bin until a signal-to-noise threshold 
is reached. After all the pixels have been accreted, it uses Voronoi 
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tessellation to make tessellated regions based on the weighted po- 
sition of the original bins. This technique has the advantage of cre- 
ating bins which are compact, varying in size smoothly with the 
surface brightness, and also provides bins with similar signal-to- 
noise ratios. We applied the method to X-ray observations of the 
Perseus cluster (Sanders et al 2004) and Abell 2199 (Johnstone et 
al 2002). Rather than use X-ray colours, we extracted spectra for 
each of the regions and used spectral fitting to derive, for example, 
temperature and abundance maps. Recently Diehl & Statler (2006) 
have generalised this algorithm to allow for data whose signal-to- 
noise does not add in quadrature. 

The motivation for further work in this area is that the methods 
above do not use the surface brightness distribution to change the 
shape of each bin. Physical parameters (e.g. density, temperature 
and abundance) usually change in the direction of surface bright- 
ness changes. The method we describe here uses the surface bright- 
ness to define bins which cover regions of similar brightness. 

Other methods have been presented for mapping the parame- 
ters of the intracluster medium. These included wavelet techniques 
(Bourdin et al 2004) and monte carlo methods (Peterson, Jernigan 
& Kahn 2004). The advantage of binning techniques is that they 
provide errors on individual spectral fit parameters, or colours, from 
a particular part of the sky. The individual measurements made us- 
ing binning methods are independent, making it easy to easily mea- 
sure the significance of individual spatial features. 

The techniques presented in this paper have already been ap- 
plied to a number of Chandra observations of clusters, including 
a deep observation of the complex structure of the Centaurus clus- 
ter (Fabian et al 2005), the possible detection of nonthermal radi- 
ation and the identification of a high metal shell likely to be asso- 
ciated with a fossil radio bubble in the Perseus cluster (Sanders 
et al 2005), a sample of moderate redshift clusters (Bauer et al 
2005), and a 900 ks observation of the Perseus cluster (Fabian et 
al 2006), finding little evidence for temperature changes associated 
with shock-like features, and producing evidence of a substantial 
reservoir of cool X-ray emitting material. 

We first present a simple smoothing method ("accumulative 
smoothing"), and then present the binning method based on the 
smoothed image ("contour binning"). 



2 ACCUMULATIVE SMOOTHING 

In order to bin using the surface brightness it was necessary for us 
to get an estimate of the surface brightness in an image in the ab- 
sence of noise and counting statistics. Simple Gaussian smoothing 
can be sufficient, but if the brightness of the object varies over its 
extent the smoothing scale will be too small or too large in parts 
of the image. There are several methods for adaptively smoothing 
an image based on the surface brightness (e.g. ASMOOTH, Ebel- 
ing et al, in preparation; CSMOOTH in CIAO, based on ASMOOTH; 
Huang & Sarazin 1996). We present a method called accumula- 
tive smoothing, which is a generalisation of the method used by 
the FTOOLS routine FADAPT, now including handling background 
and exposure variation. In addition we adapt the method to cre- 
ate smoothed colour maps. This smoothing method has the advan- 
tage of being fast, simple and easy to interpret. It allows the use 
of blank- sky background images (rather than trying to estimate the 
local background), masks, and exposure maps. Its disadvantage is 
that it does not guarantee to globally preserve counts. This is proba- 
bly not a serious problem as it is usually more important to preserve 
the local count rate. We will discuss how well it does locally later. 



It is a very simple routine which smooths an image with a top- 
hat kernel, whose size varies as a function of position. The size is 
varied so that the kernel contains a minimum signal-to-noise ratio. 

Suppose there is an input image I, containing numbers of 
counts in each pixel There is also a background image B (option- 
ally fixed to zero). The smoothed count image at positions r, where 
r is a vector representing the coordinates of the centre of the pixel 
being examined, is 
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where A' is the number of pixels summed over (~ ;ri?[r]^), and the 
smoothing radius at the pixel, R{r), is defined in Equation|2] Ej, and 
Ei are the exposure times of the foreground and background obser- 
vations, respectively. If we wish to take account of the variation in 
exposure times over an image, we can create a smoothed count rate 
image. 
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The smoothing radius R{r) is defined to be the minimum value 
of the radius (in integer units) where the signal-to-noise is greater 
or equal to a threshold value. Summing over pixels at positions a, 
where |a — r| < R{r), then the condition is met when 
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R{r) is found by incrementing R{r) from 1 until it is true. g{c) is 
an approximation for the squared-uncertainty on c counts, obtained 
from equation 7 of Gehrels (1986) when S = 1, 

g{c)=(l^f^^ . (4) 

This approximation is the larger error bar on c counts, so we over- 
estimate the negative error by assuming the errors are symmetric. 

Equation|3|is approximate if the ratio of the exposure time of 
the image and background varies over the image. We do not sum 
the squared error of B(a) as the sum of the approximation for the 
uncertainty squared g[c) will become increasingly inaccurate with 
more terms. 

Therefore the algorithm smooths the image using a top-hat 
kernel with a variable smoothing radius. The smoothing radius 
varies so that it contains a minimum signal-to-noise ratio. We han- 
dle edges and masked regions by ignoring pixels outside the valid 
region (A' does not include these pixels). 

2.1 Tests of the algorithm 

To test how closely the smoothing algorithm reproduces the sur- 
face brightness of a model image, we constructed a surface bright- 
ness model made up of six /3 surface brightness model components, 
where each component had surface brightness 

-3J3 + 1/2 



S = 5o 1 



(5) 



The components are listed in Table Q We constructed a model 
512 X 512 pixel image based on the components (Fig.Q[left]). A 
Poisson statistic realisation of the model is shown in Fig.Q(centre). 
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So (counts) 
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S/N~8 

S/N~20 

S/N~30 



Table 1. Surface brightness components. Each is a j3 surface brightness 
model (Equation |5). The core radius (r^.), component x centre (x^) and y 
centre (y^) are expressed as a fraction of the 512 x 512 pixel image. The 
coordinates are measured from the lower left of the images. 



In Fig.Q(right) we include an accumulatively smoothed image of 
the Poisson realisation. 

A good smoothing algorithm should reproduce the model im- 
age with no large scale deviations or biases, and random small 
scale deviations. We computed the fractional difference between 
smoothed images of the realisation of the model with various 
signal-to-noise threshold ratios (Fig.|2j. Also shown is a Gaussian 
smoothed fractional difference map for comparison. To aid compar- 
ison we show histograms of the fractional differences using accu- 
mulative smoothing and Gaussian smoothing with different signal- 
to-noise ratios and scale lengths in Fig.|3| 

Firstly, accumulative smoothing mostly shows deviations 
from the model which do not vary in magnitude over the image. 
The Gaussian smoothed deviations in Fig.|2|(right) increase when 
the count rate decreases. The central part of the image is uniform in 
colour, while the outer parts alternate between the extremes. Where 
accumulative smoothing shows deviations are at the edge of the im- 
age and around centrally-concentrated sources when the smoothing 
scale is too large (e.g. Fig.|2|[centre right]). Near the edges the sur- 
face brightness is too high, as there are no lower flux pixels further 
out to smooth over. If a centrally-concentrated source does not have 
enough signal to exceed the signal-to-noise threshold, M, then it is 
smoothed into its surroundings. The flux at the centre of the source 
is underestimated, and the flux in its surroundings is overestimated. 

We took a single /3 component surface brightness image (us- 
ing an image of size 1024 x 1024 pixels, and applied Gaussian 
((7 = 8 pixels using the FGAUSS program from FTOOLS) and accu- 
mulative smoothing (M = 40) to a Poisson realisation of the model. 
In the centre of the image, the Gaussian smoothing length is ap- 
proximately the accumulative smoothing length. 

The deviation distributions (Fig. |3} show that accumulative 
smoothing produces narrower tails in the distribution than Gaus- 
sian smoothing. The distributions can be well fit with a Gaussian, 
except for a tail in the positive deviations. At low signal-to-noise 
thresholds this may be due to the approximation in the uncertainties 
of the counts as symmetric. The Gehrels approximation used over- 
estimates the uncertainty on low numbers of counts, meaning that 
faint regions will have larger kernels applied than necessary when 
a low signal-to-noise threshold is used. At high signal-to-noise the 
tail is due to edge effects and oversmoothed compact sources. Us- 
ing a model distribution that contains just the central j8 component, 
and excluding the edges of the image, produces a distribution which 
is almost Gaussian (Fig.|3|[bottom]). The width of the distribution 
is 0.023, which is very close to a simple estimation of the width 
(1/M = 0.025). To get the estimate, we require c = counts in 
the absence of a background to exceed the threshold. The standard 
deviation on a mean of c counts is ^/c, or as a fraction of the mean, 
l/^=l/M. 
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Figure 3. Histogram of fractional differences between reconstructed surface 
brightness and original model. (Top) Difference between complex cluster 
model and accumulatively-smoothed data. The three curves are smoothed 
with signal-to-noise threshold ratios of 30, 20 and 8. (Centre) Fractional 
difference between the model and Gaussian-smoothed data (using widths 
of 1, 2 and 4 pixels). (Bottom) Differences between a single-component j3 
model and accumulative smoothed data (signal-to-noise threshold of 40). 
The smooth line is a Gaussian fit to the differences, with a best fitting width 
of 0.023. The edges of the image have been removed to eliminate edge 
effects. Also plotted is the result when a constant background is added to 
the model. 



© 0000 RAS, MNRAS 000, 000-000 



4 7.5. Sanders 




IS 20 50 40 20 40 60 10 20 .10 40 50 



Figure 1. (Left) Surface brightness model made up of six /3 model components. (Middle) Realisation of the model assuming Poisson statistics. (Right) 
Accumulatively smoothed image of the realisation, with a signal-to-noise ratio of 20. All three images use the same colour scale. The units are counts per 
pixel. 




Figure 2. Fractional difference between complex cluster model and accumulatively-smoothed image, for S/N 8, 20 and 30, and using Gaussian smoothing 
with a smoothing scale of 2 pixels. Note that each image uses a different colour scale to show the magnitude of the scatter. 



We have also considered the case of a simple /3 model with a 
background. We added a background of 40 counts across the entire 
image (which is the same as the peak of the j3 distribution). The 
distribution of deviations (Fig. |3| [bottom]) is fairly Gaussian, al- 
though there is a significant tail at larger deviations than the case 
without a background. In this test we assumed the background was 
very well determined (simulated by using a much larger exposure 
time for the background than the foreground). 



2.2 Colour maps 

Colour maps show the ratio of counts in two bands. Such maps are 
useful because trends in X-ray colour can follow physical trends 
(e.g. temperature or absorption). Accumulative smoothing can also 
be used to generate smoothed colour maps. Smoothing or binning 
is very useful in creating colour maps, as the error on a ratio of two 
counts can greatly vary over an X-ray image. 

Rather than considering a threshold on the signal-to-noise in 
the smoothing kernel, a maximum uncertainty on the X-ray colour 
can be used instead. 



2.3 Conclusions 



It can be concluded that accumulative smoothing, and the smooth- 
ing of FADAPT when there is no background or exposure varia- 
tion, does not introduce a major bias to an image. The median 
difference is negligible for all three smoothing thresholds shown. 
Accumulative smoothing is certainly better than Gaussian smooth- 
ing at smoothing on appropriate scales over an image. Indeed the 
deviation range does not change with position. It also produces a 
distribution of deviations that shows a narrow, almost Gaussian, 
tail, although the tail is wider for positive deviations at low signal- 
to-noise. If there is not enough signal in point sources they are 
smoothed into their surroundings, leading to deviations from the 
intrinsic surface brightness of the sky. In comparison Gaussian 
smoothing shows much more noise in regions where the count rate 
is low compared to where it is high. 



In addition the algorithm is particularly suited for presentation 
of X-ray images as an alternative to Gaussian smoothing. It may 
well be better in many cases to other forms of adaptive smoothing 
as it introduces few artifacts. 
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3 CONTOUR BINNING 

Here we describe how the accumulatively smoothed map can be 
used to define independent spatial regions. These regions can be 
used for spectral extraction and fitting, to measure physical param- 
eters, or for the calculation of X-ray colours. The method need not 
apply to accumulatively smoothed images. Other techniques such 
as adaptive smoothing could be used instead. 

Simply, the algorithm, starting at the highest flux pixel on a 
smoothed image, adds pixels nearest in surface brightness to a bin 
until the signal-to-noise ratio exceeds a threshold algorithm. A new 
bin is then created. This technique naturally creates bins which fol- 
low the surface brightness. 

3.1 Initial pass 

Contour binning takes the smoothed image, S(r). We start from the 
highest flux pixel in the smoothed image. The pixel is added to the 
set of pixels in the current bin, jS . If the signal-to-noise of the pixels 
in jS is greater or equal than a binning threshold T, 




then the region is complete. Otherwise we iterate over the subset 
of /3 which lie at the edge of the bin. We consider all the unbinned 
neighbours of these pixels and take the one which has a flux in the 
smoothed image closest to the starting flux of the bin. This pixel is 
added to /3, and we repeat the process. If there are no neighbouring 
unbinned pixels we stop this bin. We consider a pixel to be a neigh- 
bour if it differs in one of its coordinates by 1, and in the others by 
(diagonally neighbouring pixels could also be included, however). 

Once we have completed the bin, we start binning again. The 
starting pixel is again the highest flux pixel in the smoothed image 
which has not yet been binned. To optimise this potentially expen- 
sive search, we sort the pixels in the smoothed image into flux order 
before binning. 

3.2 Cleaning pass 

After all the pixels have been binned, there will be some bins with a 
signal-to-noise ratio of less than T. Many of these will contain sin- 
gle pixels which were left stranded between other bins. The num- 
ber of these stranded bins decreases with increased smoothing of 
the initial image. We therefore "clean up" the bins to remove any 
below the threshold. We start by selecting the bin with the lowest 
signal-to-noise ratio. By examining its edge pixels, we select the 
pixel which is closest in smoothed value to a pixel in a neighbour- 
ing bin. This pixel is moved into the corresponding neighbouring 
bin (thereby increasing the signal-to-noise ratio of its neighbour). 
We repeat, transferring the pixels until there are no pixels remaining 
in the lowest signal-to-noise bin. We move onto the next remaining 
bin which has the lowest signal-to-noise and repeat the whole pro- 
cess. The cleaning pass ends when there are no more bins with a 
signal-to-noise ratio of less than the threshold T . 

3.3 Constraints 

On a high signal-to-noise dataset with low signal-to-noise bins the 
above process is sufficient. However there are possible undesirable 
features under other conditions. If the smoothed map is not very 



smooth, the bins can become distended (this can be improved at 
the detriment of potential detail by increasing M). If the smoothed 
map is very smooth and radial, then bins will become whole annuli 
if T is too large. In some cases this outcome is desirable. 

To prevent these problems we can impose extra geometric 
constraints during the binning process. Ideally such a constraint 
could be computed in constant time on each iteration or else grow- 
ing large bins becomes prohibitively expensive. We identified a fast 
constraint that gave the most aesthetic results. The procedure is to 
estimate a radius R that the current bin would have if it were a cir- 
cle, i.e. R --JN 1%, where A' is the number of pixels in the bin. 
By examining circles on a pixelated grid, we calculate this value 
exactly for integer radii. This is far better than the simple esti- 
mate for small circles. We then measure the distance d of the pixel 
to be added from the current flux-weighted centroid of the bin. If 
Rjd > C, where C is the constraint parameter, then we do not add 
the pixel. The constraint is also applied during the cleaning pass 
by moving pixels into only bins in which the constraint would be 
satisfied. If there are no bins which would satisfy the constraint, the 
constraint is broken to prevent isolated pixels. 

3.4 Tests of the algorithm 

We binned the simulated image in Fig.Qwith the algorithm. Fig.|4| 
shows the image binned using a signal-to-noise ratio, T threshold of 
30 and 100 (smoothing with ratios, M, of 15 and 40, respectively). 
The r = 100 image has also had a constraint of C = 2 applied. We 
also show a binned image using a signal-to-noise ratio of 100 us- 
ing the bin-accretion algorithm of Cappellari & Copin (2003), and 
simple adaptive binning (Sanders & Fabian 2002) for comparison. 
The contour binning method appears to give a good reconstruction 
of the surface brightness. Using the same signal-to-noise ratio, it 
follows the surface brightness contours better than the tessellated 
image. 

In Fig.|5|we show absolute fractional difference maps between 
the binned images in Fig. |4|and the model map (Fig. Q left). The 
comparison between the contour binning and bin accretion meth- 
ods is not completely fair, as the contour binning procedure uses a 
minimum signal-to-noise ratio of 100, whilst bin accretion uses a 
mean signal-to-noise ratio of 100. However increasing the signal- 
to-noise ratio used by the bin accretion method would increase the 
size of the bins further. 

The signal-to-noise ratio of the bin accretion method can be 
adjusted to give almost the same number of bins as contour bin- 
ning. The distribution of the deviations of pixels from the model 
are shown for the three methods in Fig.|6| The plot shows that con- 
tour binning produces values which are considerably closer to the 
model distribution than the other methods, due to its method of 
grouping together pixels with the similar values. 

For contour binning, the deviations appear random over the 
bins, except for the two peaked sources with large signal-to-noise 
ratios. In Fig. |7|is shown the distribution of the fractional differ- 
ences between the binned surface brightness map and the model 
map for the contour binned images shown, plus another with T = 
10. The plots show the the differences are close to Gaussian. At 
the low signal-to-noise ratio of 10 there is some asymmetry in the 
distribution towards higher fractional differences. This may be due 
to the use of the Gehrels approximation for the Poisson errors. The 
width of the distributions is close to the estimated value of 1 /T. The 
estimated width of the distribution is the same as at the end of Sec- 
tion l2.1l c = counts are required to exceed the threshold without 
a background. The standard deviation on a mean of c counts is ^/c. 
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Figure 4. (Left) Contour binned surface brightness image with T = 30, M = 15. (Centre Left) Contour binned image with T = 100, A/ = 40 and constraint 
C = 2 applied. (Centre Right) Image created using bin-accretion algorithm of Cappellari & Copin (2003) with signal-to-noise ratio of 100. (Right) Adaptively 
binned image (Sanders & Fabian 2001) with a fractional error of 0.01 (equivalent to a signal-to-noise ratio of 100). All images have units of counts per pixel 
and the colour scales are the same as in Fig.fTI 






Figure 5. Absolute fractional differences between binned surface brightness image and model image. The panels shown are the same binning as Fig.l4l 



If this is turned into a fractional width of the distribution it equals 

i/7?=i/r. 



3.5 Spectral fitting tests 

To thoroughly test the algorithm when using it to define regions to 
perform spectral extraction and fitting on, we constructed simulated 
Chandra event files for a cluster of galaxies. To produce the files, 
we iterated over cells in a 3-dimensional cube of dimensions 500 x 
500 X 1000 pixels, using a cell size in {x.,y,z) of 2, 2 and 4 pixels. 
In each cell we used XSPEC to fake a PI (pulse-invariant) spectrum 
corresponding to the gas in the cluster that would be in that cell. 
To do this we use an analytical form of the temperature, density 
and abundance of the gas as a function of radius. The density is 
used to calculate the emission measure, and the spectrum is faked 
using the MEKAL spectral model (Mewe, Gronenschild & van den 
Oord 1985; Liedahl, Osterheld & Goldstein 1995), with Galactic 
photoelectric absorption. 

We added the spectra along each line of sight in z to produce a 
2-dimensional set of spectra. For each spectrum, we manufactured 
faked events which would produce the spectrum if it were extracted 
from the 2x2 pixel projected cell on the sky. The events were 
randomised in position over the projected cell. 

In detail, a template event file from a real Chandra observa- 
tion was populated with these events. To create the faked spectra 
a constant response and ancillary-response was used in each cell. 




0.5 0.75 

Fractional deviation from model image 

Figure 6. Plot of the distribution of absolute deviation of binned pixels 
from the model values for similar signal-to-noise ratios. The contour bin- 
ning and adaptive binning values were taken from Fig.l5l(centre left and far 
light panels). The signal-to-noise ratio of the bin accretion used has been 
adjusted from 100 to 106 to get a similar number of bins compared to the 
contour binning results (257 versus 255). This makes little difference to the 
distribution. The plot shows that contour binned surface brightness is closer 
to the model surface brightness than the other methods. 
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Figure 7. Distribution of fractional difference between the binned surface 
brightness of a bin and the mean model suii'ace brightness in the bin. The 
smooth curve is a best fitting Gaussian. (Top) Fractional differences for 
T = 100, M = 40, constraint applied with C = 2. Best-fitting Gaussian cen- 
tre is 3.0 X lO""*, width is 9.4 x 10"^ (Centre) T = 30, A/ = 15. Gaussian 
parameters of (-6.8 x 10"", 0.037) (Bottom) r = 10, M = 15. Gaussian 
parameters of (1.3 x 10"^ 0.087). 



We used a response for the centre of the ACIS-S3 CCD in 2000. 
The spectra were faked for an observation of 60 ks. In addition to 
the emission of the cluster in each cell on the sky, we added events 
which would produce a background spectrum corresponding to the 
spectrum observed from a blank-sky observation (generated by fit- 
ting a three-component power law model to spectra from a blank- 
sky background field). 

To model the emission from a cluster, we fitted analytic mod- 
els to the deprojected temperature, density and abundance profiles 
measured in Abell 2204 (Sanders, Fabian and Taylor 2005a). The 
cluster was centred in the centre of the cube, and the properties of 
the gas at each radius in the cluster was found using the analytic 
fits to the A2204 profiles. 

The density profile in A2204 was fitted with a j8 density model 
of functional form 



-3P/2 



with r,. = 13.97 kpc, no = 0.1978 cm- 
perature was modelled with 



r(r) = ro+ri 



(7) 



and P = 0.465. The tem- 



(8) 



with To = 2.98 keV, Ti = 7.38 keV, = 50.2 kpc and rj = 2.97. 
The abundance was modelled with the same functional form, with 
Zo = O.156Z0, Zi = 0.921, rc = 47.8 kpc, and 77 = -7.56. The 
simulated cluster was assumed to lie at a redshift of 0.15, and we 
used a cosmology of Ho = 70 km s" ' Mpc" ' and £1/^ = 0.7. Using 
Chandra pixel sizes, 1 pixel corresponds to 1.28 kpc for the real 
cluster. 

Using the faked dataset, we created an image of the simu- 
lated cluster, and of the simulated background, in the 0.5 to 7.0 
keV band. This was contour binned using a signal-to-noise of 75 
(--^ 5600 counts), smoothing the image with a signal-to-noise of 20, 
and constraining the bins using a parameter C = 2 (see Section l3!3l . 
This procedure yielded 52 bins in total. 

Individual spectra were extracted from the dataset using the 
regions defined by each bin, and each spectrum was fitted with a 
MEKAL model with free absorption, metallicity, temperature and 
normalisation. The spectra were fit between 0.5 and 7 keV, min- 
imising the statistic. We show the resulting maps in Fig. |8| 
whilst radial profiles are shown in Fig.|9| 

The results show the accumulative binning algorithm does a 
reasonable job at reproducing the input profiles. They do demon- 
strate some deviations which are due to projection effects, and 
emission weighting, however. Any method which does not account 
for projection will be susceptible to these problems. In particular 
the cool gas in the centre is measured to be too hot. This is due to 
the projection of the hotter gas outside this region. To account for 
projection, you must assume some sort of symmetry. 

As a comparison, Fig. llOl shows the temperature map from the 
data binned using the bin accretion algorithm, with the same signal- 
to-noise ratio. In the central region the bins are not as well shaped to 
match the surface brightness, leading to a blocky appearance. There 
are some problems with small bins being generated in the outskirts 
of the image, but it is possible this could be due to problems in our 
implementation of the algorithm. Again the comparison is not quite 
fair as bin accretion uses a mean signal-to-nose ratio, rather than a 
minimum one. However, increasing the signal-to-noise ratio would 
make the bins larger. 
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Figure 8. Resulting emission weighted temperature, pliotoelectric absorption and abundance maps generated from the simulated cluster Each region was 
created to have a minimum signal-to-noise ratio of 75. Fig.|9]shows the extracted radial profiles. 
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Figure 9. Profiles produced by fitting bins with a minimum signal-to-noise ratio of 75 using the faked data. (Left) temperature, (centre) photoelectric absorption 
and (right) abundance. The solid lines show the intrinsic (not projected) profiles the faked data was generated with. 



4 CASSIOPEIA A 

Cassiopeia A is a very good target to demonstrate the algorithm 
on. The large range in brightness of the image between the bright 
knots and the darker regions shows off the advantages of this par- 
ticular method. We present here results using a 49.7 ks observation 
of Cas A using the ACIS-S3 detector on Chandra taken on 2002- 
02-06 (observation number 1952). 



4.1 Data preparation 

Firstly we reprocessed the level 1 event file using CIAO 3.3.0.1 with 
the gain file acisD2000-08-12gainN0005.fits. Time dependent gain 
correction was also applied using the CORR.TGAIN utility using the 
November 2003 correction (Vikhlinin 2003). In addition the CIAO 
version of the LC_CLEAN script (Markevitch 2004) was used to re- 
move periods of the dataset where the count rate was 20 per cent 
away from the mean. Very little time was removed, producing a 



level 2 event file with an exposure of 48.3 ks, although such filter- 
ing is difficult on such a bright target. 

The dataset was taken using the unusual GRADED observa- 
tion mode of the ACIS detector, where events likely to be back- 
ground are discarded at the satellite rather than removed by the 
observer. Despite this, the background rate at high energy (9 keV 
or more) matches that from blank sky observations. We used a tai- 
lored version of a blank sky observation to account for background 
in this analysis. 

As Cas A is a very bright object it is important to take account 
of 'out of time' events. These are events which occur whilst the de- 
tector is being read out. As the ACIS detectors use framebuffers this 
time is short. Around 1.3 per cent of events occur out of time. We 
used the MAKE_READOUT_BG script (Markevitch 2003) to generate 
a synthetic event file to include in the background subtraction. The 
script works by randomising the CHIPY coordinate of the events in 
the input event file (which is the actual observation), and increasing 
the exposure time by the reciprocal of out of time fraction. In this 
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Figure 10. Temperature map created by applying bin accretion algorithm 
to fake data. 
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Figure 13. Histogram of fractional deviations of each pixel of a binned im- 
age from accumulatively smoothed image of Cas A (signal to noise thresh- 
old of 8). The area examined was limited to a 2.8 arcmin radius around the 
core of the remnant. Shown are the results for contour binning, adaptive 
binning and bin accretion using similar signal-to-noise ratios. 
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Figure 11. Background subtracted, accumulatively smoothed image of 
Cas A in the 0.5 to 7 keV band, smoothed using a signal-to-noise ratio 
threshold, M of 15. 



observation the out of time background dominates over the blank 
sky background over much of the energy range. 



4.2 X-ray image 

In Fig. ^3 is shown an accumulatively smoothed X-ray image of 
the dataset. We show in Fig. 1121 a contour binned image of the 
dataset with a signal-to-noise threshold, T, of 100 (the image was 
smoothed with a signal-to-noise threshold of 15 before binning) 
along with adaptively and accretion binned images with similar 
signal-to-noise ratios. The binning procedure created 1176 spatial 



bins, giving an average of 1.2 x lO'* counts per spectrum before 
background subtraction. It can be seen that the bins follow the 
shape of the surface brightness well, and better than the other bin- 
ning methods. Note that there is relatively bright emission around 
the remnant (this is difficult to see in plots shown, due to the 
greyscale colour scaling). This emission is in fact radiation from 
Cas A scattered by dust in the interstellar medium. The measure- 
ments presented later include results outside of the remnant which 
are from the scattered radiation. 

To compare the binning methods better, we can compare the 
surface brightness of binned pixels against what they are expected 
to be. If we take the absolute fractional difference between a binned 
accumulately smoothed image of the remnant and an accumula- 
tively smoothed image, then plot the histogram of the differences, 
we obtain the distributions in Fig. ll3l The area examined was lim- 
ited to a radius of 2.8 arcmin from the centre of the remnant to avoid 
edge effects. The plot shows that contour binning is much better at 
reproducing the accumulately smoothed surface brightness of the 
object than the other methods. This argument assumes that an ac- 
cumulative smoothed image is a good estimate of the intrinsic sky 
surface brightness. To make the comparison fairer we increased the 
signal-to-noise used by the bin accretion method to 109 to match 
the number of bins used by contour binning, although it makes very 
little difference to the results. The number of bins used by contour 
binning, bin accretion and adaptive binning are 1 175, 1 182 and 889 
respectively. Contour binning also produces the most accurate sur- 
face brightnesses from our simulated data (Fig. |6|. The test with 
the simulated data is less biased as we know the intrinsic surface 
brightness of the object, although the results are similar. 

4.3 Spectral fitting 

We extracted spectra from the events file in each of the binned 
regions. Rather than use the standard CIAO DMEXTRACT tool we 
used our own code which extracts all the spectra in a single pass. 
The spectra were grouped to have a minimum number of 20 counts 
in each spectral channel using GRPPHA from FTOOLS. Using the 
CIAO MKWARF and MKRMF tools, we generated weighted response 
and ancillary response files for each of the regions, weighted by the 
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Figure 12. Binned images of Cas A. (Left) Contour binned image using a signal-to-noise binning thi'eshold ratio, T, of 100. The image was smoothed using 
M = 15 before binning. (Centre) Adaptively binned image with a fractional error threshold of 0.01 using algorithm of Sanders & Fabian (2001). Note that this 
binning did not take account of the readout background. (Right) Bin accretion binned image with a signal-to-noise of 100 using algorithm of Cappellari & 
Copin (2003). 



number of counts between 0.5 and 7 keV. In addition we extracted 
background spectra from the blank sky and out of time event files. 

The spectra were fit using XSPEC 11.3.2 between energies of 
0.6 and 7 keV. We fitted a VNEI model to each region, with the 
neivers parameter set to 2.0, which uses the ionisation fractions of 
Mazzotta et al (1988) and uses APED (Smith et al 2001) to gen- 
erate the resulting spectrum. The emission model was absorbed 
by a PHABS absorption model (Balucinska-Church & McCammon 
1992). In the fit the absorption, temperature, normalisation, ioni- 
sation timescale, redshift, Ne, Mg, Si, S, Ar, Ca, Fe and Ni abun- 
dances were free parameters (using a maximum of 40 Zp, in each 
abundance parameter). The H, He, C, N and O abundances were 
fixed at solar values. We used the solar abundance ratios of Anders 
& Grevesse (1989). In some regions there are prominent lines at 
Fe-K energies (e.g. Hwang et al 2000), which lead to poor reduced 
values. We therefore added a Gaussian component to the model 
allowed to vary between 6 and 7 keV. 

With the large number of free parameters it was easy for the 
fits to become stuck in local minima. We tried to avoid this by free- 
ing the most influential parameters first, and doing a fit after each 
freeing. Also during the fit we searched for errors on the parameters 
to help find other minima. However, this prescription does not 
guarantee to find the absolute minimum in x^ space, so some pa- 
rameters may not be at their best-fitting values. In addition the best 
fitting parameters may not be physically sensible. After the spec- 
tra were fit, the best-fitting values of each parameter were used to 
generate maps. 

The model used is somewhat simplistic compared to those 
used by other authors, not taking account of multi-temperature re- 
gions or allowing the lines from each element to have different ve- 
locities (except for Fe-K). Most of the features reported by oth- 
ers are reproduced, however, as we discuss below. The model also 
assumes the emission is thermal in origin, and does not include a 
non-thermal component. Hughes et al (2000) have noted continuum 
dominated regions. These regions can be fitted by a powerlaw com- 
ponent or an absorbed bremsstrahlung. Hughes et al (2000) prefer 
the interpretation that the continuum is X-ray synchrotron emis- 
sion. 



In Fie. 1141 are shown the best fitting temperature, absorbing 
column density, ionisation timescale and velocity. The velocity was 
obtained by multiplying the best fitting redshift by the speed of 
light. Abundances for Si-like elements are shown in Fig. 1151 and for 
Fe-like elements in Fig.^| We finally display in Fig. ll7l the nor- 
malisation of the main VNEI component per unit area, the reduced- 
X^ of the fits, the normalisation of the Gaussian Fe-K component 
and its energy. All the maps here have been smoothed with a Gaus- 
sian of 2 pixels (0.98 arcsec) for display purposes. The positions in 
the map and text are in J2000 coordinates. 

4.4 Discussion 

The absorption map (Fig. 1 141 [top right]) shows a dense clump of 
X-ray absorbing material to the west of the nebula in great detail. 
This material was itself found by previous observations (Keohane, 
Rudnick & Anderson 1996; Willingale et al 2002). 

The temperature map shows the location of the inner and outer 
shocks well (Gotthelf et al 2001), as regions of high temperature, in 
addition to the direction of the north-east jet. The spectral fits also 
indication the position of the south-west jet (Hwang et al 2004) as 
higher temperature regions. 

The abundance maps can be compared against those created 
by plotting the ratio of emission in a band crafted to the emis- 
sion line of a particular element to continuum emission (Hwang, 
Holt & Petre 2000). This technique can only be used for elements 
which produce lines easily identifiable from the spectrum, however. 
Our maps show good agreement with the morphology of the ratio 
maps, although are some point-to-point differences. For example, 
we find more emission from S and Ar where the northern enhance- 
ment in abundance meets the eastern region, around (23*'33™33'* 
+48° 50' 05"). In addition, although the Hwang et al (2000) maps 
show a rim of Si (and possible S) enhancement to the south of the 
object, we see the rim break up into separate peaks (as found in 
the megasecond observation of Hwang et al 2004), and see the en- 
hancement in S, Ar, Ca, less strongly in Fe, and possibly in Mg and 
Ne. Again Si, S, Ar and Ca vary much more uniformly than Ne, 
Mg, Fe and Ni. 
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Although we do not allow separate velocities for the differ- 
ent metal components (except Fe-K), the velocity map (Fig. 1141 
bottom right), where velocities towards the observer are shown as 
negative, shows similar features to those found by others for the Si 
line which is usually dominant (e.g. Willingale et al 2002). We are 
able to map the velocity in those regions where the line emission is 
weak, as the emission is binned over larger regions. The absolute 
velocity normalisation appears to be uncertain, as earlier calibration 
versions with the same data gave velocities systematically around 
3500 km s^' larger. Using the Fe-K line energy (Fig. 1171 bottom 
right), the trend of velocity is again very similar to the Willingale 
et al. results, but here there is not enough signal to map the energy 
in low intensity regions. 

The plot of the normalisation of the VNEI component is par- 
ticularly striking (Fig. ll7l top left). It shows the emission from the 
remnant when the line emission has been removed. Cas A is much 
more symmetrical in the underlying emission. 

The reduced-;!^^ is not too bad generally when the Fe-K emis- 
sion has been modelled (as it has here), given the complexity of the 
spectra and the simplicity of the fitted model. Where the reduced- 
is at its worst, the model seems unable to take account of the 
continuum and the emission lines together. Adding a second tem- 
perature component in these regions improves the quality of the fit 



substantially. Multiple temperature fits have been necessary before 
in the spectral fitting (e.g. Willingale et al 2002). There is much 
room for further investigation of the spectra here, but it is not the 
aim of this paper to study the physics of these regions in detail. 



5 CONCLUSIONS 

We have shown that the accumulative smoothing method, or that 
used by FADAPT, is a robust way to estimate the real surface 
brightness distribution of an extended object. We have discussed 
the method in the cases of varying exposures, and including back- 
grounds. 

We have described a new method for choosing regions for 
spectral extraction, or colour analysis, called contour binning. We 
have shown that method reliably creates bins which follow the sur- 
face brightness. We demonstrate the method using a simulated clus- 
ter dataset, and apply it to a Chandra observation of the supernova 
remnant Cassiopeia A. 

The algorithm is ideal where spectral changes are associated 
with changes in surface brightness, as is often the case in X-ray 
astronomy. The results on objects with fine spatial detail are very 
closely matched to the object, aesthetically pleasing, and easy to 
interpret, compared to other methods (as exemplified by Fig. ll2l and 
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Figure 15. Abundances in Cas A of Si, S, Ar and Ca relative to solar. 



Fig. ll3> . As the algorithm follows the surface brightness it naturally 
bins regions which are likely to be physically associated. It avoids 
mixing spectra together for physically disassociated regions. 



AVAILABILITY 

A C++ implementation of the algorithm can be downloaded 
from http: / /www-xray . ast . cam. ac.uk/papers/contbin/J , 
including other helpful associated programs. 
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Figure 16. Abundances in Cas A of Ne, Mg, Fe and Ni relative to solar. 
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Figure 17. VNEI normalisation per unit area (a proxy for the emission when the hne emission has been removed), reduced of the fit, normahsation of the 
Fe-K line per unit ai'ea, and energy of the Fe-K line in Cas A. 
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